#####################################
#### Assign necessary libraries. ####

starttime <- Sys.time()

library(survival)
library(plyr)
library(optmatch)
library(rcbalance)
library(rcbsubset)
library(car)
library(data.table)
library(splitstackshape)
library(RCurl)
library(MASS)
library(pROC)


source("/Users/sharpe/Desktop/Young Surgeons/code/apply_penalty.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/subtabfun.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/meantab_both.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/baltab_ys.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/Source/smahal_func.R")
source("/Users/sharpe/Desktop/Young Surgeons/code/Source/penalize.near.exact.v2.R")



#####################################################
#### Read in all of the datasets. ###################

trad <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_trad_ortho_simul_11_02.csv")
mod <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_mod_ortho_simul_11_02.csv")



################################
#### Important Information. ####
#### newID_1 = BENE_ID.
#### newID_2 = AT_PHYSN_NPI.
#### newID_3 = AT_PHYSN_UPIN.
#### newID_4 = OP_PHYSN_NPI.
#### newID_5 = OP_PHYSN_UPIN.
#### newID_6 = PRVDR_NUM.
#### newID_7 = UPIN.
#### newID_8 = NPI.
#### newID_9 = ResearchID.
#### newID_10 = CLM_ID.


##########################################
#### Set a seed for random processes. ####

set.seed(11292017)


###########################################################################################################################################
#### Prepare the narrow procedure groups. #################################################################################################
#### Create the spx_none variable. ########################################################################################################
#### Create a flag for whether a patient was seen before Oct. 2010. #######################################################################
#### Create binary admit year variables. ##################################################################################################
#### Create sum of comorbidities variable. ################################################################################################
#### Create interaction variables in order to penalize only the traditional or modern era patients later on for certain comorbidities. ####


comorbidities <- c("silb_CHF","silb_PastMi","silb_PastARR","silb_ANGINA","silb_Valvulardis","silb_UnstAngina","silb_HTN","silb_DEMENT","silb_LIVER","silb_PARAPLEG","silb_RENLDYS","silb_renlfail",
                   "silb_SEIZUR","silb_STROKE","silb_AbdCancer","silb_ASTHMA","silb_CANCER","silb_CHRNLUNG","silb_DIABETES","silb_PostPulmFibr",
                   "silb_AIDS","silb_ChrPepticUlcer","silb_COAG","silb_CollagenVascular","silb_CongenitalCoag","silb_CUSHINGS","silb_GRAVES","silb_HYPOTHY","silb_SickleCell","silb_THROMB")

prep_npgs <- function(data){
  
  data$Narrow_Procedure_Group <- gsub(pattern=" ", replacement="_", x=data$Narrow_Procedure_Group) 
  data$Narrow_Procedure_Group <- gsub(pattern="-", replacement="", x=data$Narrow_Procedure_Group)
  data$Narrow_Procedure_Group <- gsub(pattern="\U3e36393c_", replacement="", x=data$Narrow_Procedure_Group)
  data$Narrow_Procedure_Group <- gsub(pattern="\\(", replacement="", x=data$Narrow_Procedure_Group)
  data$Narrow_Procedure_Group <- gsub(pattern="\\)", replacement="", x=data$Narrow_Procedure_Group)
  data$not_oct_2010_plus <- ifelse(data$oct_2010_plus==1,0,1)
  
  data$year_1999 <- ifelse(data$admit_year==1999,1,0)
  data$year_2000 <- ifelse(data$admit_year==2000,1,0)
  data$year_2001 <- ifelse(data$admit_year==2001,1,0)
  data$year_2002 <- ifelse(data$admit_year==2002,1,0)
  data$year_2003 <- ifelse(data$admit_year==2003,1,0)
  data$year_2009 <- ifelse(data$admit_year==2009,1,0)
  data$year_2010 <- ifelse(data$admit_year==2010,1,0)
  data$year_2011 <- ifelse(data$admit_year==2011,1,0)
  data$year_2012 <- ifelse(data$admit_year==2012,1,0)
  data$year_2013 <- ifelse(data$admit_year==2013,1,0)
  
  data$parapleg_mod <- data$era_mod*data$silb_PARAPLEG
  data$congcoag_mod <- data$era_mod*data$silb_CongenitalCoag
  data$collvasc_mod <- data$era_mod*data$silb_CollagenVascular
  data$collvasc_trad <- data$era_trad*data$silb_CollagenVascular
  data$hypo_trad <- data$era_trad*data$silb_HYPOTHY
  
  data$comorb_sum <-  apply(data[,c(comorbidities)], MARGIN=1, FUN=sum)
  data$renal_fail_trad <- data$era_trad*data$silb_renlfail
  
  data$past_adm_trad <- data$past_adm_6mo*data$era_trad
  data$stroke_trad <- data$silb_STROKE*data$era_trad
  data$parapleg_mod <- data$silb_PARAPLEG*data$era_mod
  data$seizure_mod <- data$silb_SEIZUR*data$era_mod
  data$unstang_mod <- data$silb_UnstAngina*data$era_mod
  data$hypo_mod <- data$silb_HYPOTHY*data$era_mod
  
  
  return(data)
  
}


trad <- prep_npgs(data=trad)
mod <- prep_npgs(data=mod)





#################################################
#### Get how many patients each surgeon saw. ####

all_patients <- rbind(trad,mod)

n_patients_all <- as.data.frame(table(all_patients$newID_9))
colnames(n_patients_all) <- c("newID_9", "n_patients")

trad <- merge(x=trad, y=n_patients_all,by="newID_9",all.x=TRUE)
mod <- merge(x=mod, y=n_patients_all, by="newID_9", all.x=TRUE)


####################################
#### Create a new surgeon flag. ####

trad$new_surgeon_flag <- ifelse(trad$Trad_New_Surgeon_flag_numeric==1,1,0)
mod$new_surgeon_flag <- ifelse(mod$Mod_New_Surgeon_flag_numeric==1,1,0)


####################################################
#### Get a list of all narrow procedure groups. ####

names_tg <- as.character(unique(trad$Narrow_Procedure_Group))
names_mg <- as.character(unique(mod$Narrow_Procedure_Group))

all_names <- c(names_tg, names_mg)
all_names <- unique(all_names)
all_names_char <- as.character(all_names)

nl <- length(all_names)

create_binary_npgs <- function(data){
  
  start_n_vars <- dim(data)[2]
  
  for(i in 1:nl){
    
    vn <- paste(all_names[i])
    data$temp_var <- ifelse(data$Narrow_Procedure_Group==all_names[i], 1, 0)
    current <- i+start_n_vars
    colnames(data)[current] <- vn
    
  }
  
  
  start_n_vars <- dim(data)[2]

  for(i in 1:nl){

    vn <- paste(all_names[i])
    data$temp_var <- data[,vn]*data[,"oct_2010_plus"]
    current <- i+start_n_vars
    colnames(data)[current] <- paste(vn, "2010_plus",sep="_")


  }

  start_n_vars <- dim(data)[2]

  for(i in 1:nl){

    vn <- paste(all_names[i])
    data$temp_var <- data[,vn]*data[,"not_oct_2010_plus"]
    current <- i+start_n_vars
    colnames(data)[current] <- paste(vn, "before_2010",sep="_")


  }

  return(data)
  
}

trad <- create_binary_npgs(data=trad)
mod <- create_binary_npgs(data=mod)


#######################################################
#### Randomly sample 10 of the new surgeons cases. ####

trad_new <- subset(trad, Trad_New_Surgeon_flag_numeric==1)
mod_new <- subset(mod, Mod_New_Surgeon_flag_numeric==1)


trad_exp <- subset(trad, Trad_Exp_Surgeon_flag_numeric==1)
mod_exp <- subset(mod, Mod_Exp_Surgeon_flag_numeric==1)

exp_patients_before <- rbind(trad_exp, mod_exp)
exp_patients_before$treatment <- 0


n_pats_exp <- as.data.frame(table(exp_patients_before$newID_9))
colnames(n_pats_exp) <- c("newID_9","n_patients")

exp_patients_before <- merge(x=exp_patients_before, y=n_pats_exp, by="newID_9")

NEW_PATIENTS_BEFORE_TRAD <- trad_new
NEW_PATIENTS_BEFORE_MOD <- mod_new

NEW_PATIENTS_BEFORE_TRAD$treatment <- 1
NEW_PATIENTS_BEFORE_MOD$treatment <- 1

#### Save the patients that were not randomly sampled to be used in the risk score. #######################################################################################################
#### sampled_data_trad and _mod are lists, where the first level is a dataset of the selected 10 patients per new surgeon and the second level are the patients who were not selected. ####
sampled_data_trad <- stratified(indt=trad_new, group=c("newID_7"), size=10, bothSets=TRUE)
sampled_data_mod <- stratified(indt=mod_new, group=c("newID_8"), size=10, bothSets=TRUE)

trad_match_patients <- as.data.frame(sampled_data_trad[[1]])
trad_risk_patients <- as.data.frame(sampled_data_trad[[2]])

mod_match_patients <- as.data.frame(sampled_data_mod[[1]])
mod_risk_patients <- as.data.frame(sampled_data_mod[[2]])

trad <- rbind(trad_match_patients, trad_exp)
mod <- rbind(mod_match_patients, mod_exp)

NEW_PATIENTS_BEFORE <- rbind(trad_match_patients, mod_match_patients)

n_pats <- as.data.frame(table(NEW_PATIENTS_BEFORE$newID_9))
colnames(n_pats) <- c("newID_9","n_patients")

NEW_PATIENTS_BEFORE <- merge(x=NEW_PATIENTS_BEFORE, y=n_pats, by="newID_9")



#########################################################################################################################################
#### Create a propensity score for the randomly sampled patients from new surgeons, and all of the patients of experienced surgeons. ####
#### Do this separatetly for each era. ##################################################################################################
#### UPDATE: We are no longer using a propensity score at the surgeon level. ############################################################

all_names_no_cpl <- all_names[which(!(all_names %in% c("colectomy_partial__lap")))]
hip_npgs <- all_names[grep(pattern="hip", x=all_names)]
knee_npgs <- all_names[grep(pattern="knee", x=all_names)]
ortho_npgs <- c(hip_npgs, knee_npgs)
ortho_npgs_for_prop_score <- ortho_npgs[-length(ortho_npgs)]

gensurg_npgs <- all_names_no_cpl[which(!(all_names_no_cpl %in% ortho_npgs))]
gensurg_npgs_for_prop_score <- gensurg_npgs[-length(gensurg_npgs)]


###################################################################################
#### Create a dataset of the experienced surgeons and one of the new surgeons. ####
#### Create a flag for treatment group (needed for the match). ####################


trad_new_surgeons <- subset(trad, Trad_New_Surgeon_flag_numeric==1)
trad_exp_surgeons <- subset(trad, Trad_Exp_Surgeon_flag_numeric==1)


mod_new_surgeons <- subset(mod, Mod_New_Surgeon_flag_numeric==1)
mod_exp_surgeons <- subset(mod, Mod_Exp_Surgeon_flag_numeric==1)


trad_new_surgeons$treatment <- 1
mod_new_surgeons$treatment <- 1

trad_exp_surgeons$treatment <- 0
mod_exp_surgeons$treatment <- 0

trad_new_exp_surg <- rbind(trad_new_surgeons, trad_exp_surgeons)

mod_new_exp_surg <- rbind(mod_new_surgeons, mod_exp_surgeons)



pre_surgeon_match <- function(data){
  
  surgeon_npg_sums <- aggregate(x=data[,c(ortho_npgs, paste(ortho_npgs, "before_2010",sep="_"), paste(ortho_npgs, "2010_plus",sep="_"))],by=list(data$newID_9),FUN=sum)
  names(surgeon_npg_sums)[1] <- "newID_9"
  names(surgeon_npg_sums)[2:length(names(surgeon_npg_sums))] <- paste("sum", names(surgeon_npg_sums)[2:length(names(surgeon_npg_sums))],sep="_")
  assign("surgeon_npg_sums", names(surgeon_npg_sums)[2:length(names(surgeon_npg_sums))], envir=.GlobalEnv)
  
  data <- merge(x=data, y=surgeon_npg_sums, by="newID_9", all.x=T)
  
  
  oct_sums <- aggregate(x=data[,c("oct_2010_plus","not_oct_2010_plus")], by=list(data$newID_9), FUN=sum)
  names(oct_sums)[1] <- "newID_9"
  names(oct_sums)[2:length(names(oct_sums))] <- c("oct_2010_plus_sum","not_oct_2010_plus_sum")
  data <- merge(x=data, y=oct_sums, by="newID_9", all.x=T)
  
  #### For each surgeon, calculate the mean exp_at_partb. ####
  
  mean_exp_at_partb <- aggregate(x=data$exp_at_partb, by=list(data$newID_9), FUN=mean)
  colnames(mean_exp_at_partb) <- c("newID_9", "mean_exp_at_partb")
  
  data <- merge(x=data, y=mean_exp_at_partb, by="newID_9", all.x=T)
  
}


trad_new_exp_surg <- pre_surgeon_match(data=trad_new_exp_surg)
mod_new_exp_surg <- pre_surgeon_match(data=mod_new_exp_surg)

#########################################################
#### Create datasets that contain 1 row per surgeon. ####

unique_surgeons_trad <- trad_new_exp_surg[!duplicated(trad_new_exp_surg$newID_9),]
unique_surgeons_mod <- mod_new_exp_surg[!duplicated(mod_new_exp_surg$newID_9),]


length(unique(trad_new_exp_surg$newID_9))

#########################################
#### Get number of unique hospitals. ####

trad_new_exp_surg <- trad_new_exp_surg[order(-trad_new_exp_surg$treatment),]
hospitals_trad <- unique(trad_new_exp_surg$newID_6)
nst <- subset(trad_new_exp_surg, treatment==1)
est <- subset(trad_new_exp_surg, treatment==0)

mod_new_exp_surg <- mod_new_exp_surg[order(-mod_new_exp_surg$treatment),]
hospitals_mod <- unique(mod_new_exp_surg$newID_6)
nsm <- subset(mod_new_exp_surg, treatment==1)
esm <- subset(mod_new_exp_surg, treatment==0)

############################################################################
#### Create an empty list to store the Mahalanobis structure. ##############
#### SWITCH TO USING ResearchID (newID_9), so we can follow surgeons. ######

surgeon_level_match_trad <- vector(mode="list", length=length(unique(est$newID_9)))
surgeon_level_match_mod <- vector(mode="list", length=length(unique(esm$newID_9)))

#####################################################################
#### Perform the new-to-experienced surgeon match for both eras. ####

source("/Users/sharpe/Desktop/Young Surgeons/code/Source/new_exp_surg_match_counts_w_exact_oct_and_proc.R")
trad_match <- NULL
mod_match <- NULL

trad_match <- new_exp_surg_match_counts(data=as.data.frame(unique_surgeons_trad), hospid="newID_6", hospitals=hospitals_trad, surg_id_q ="newID_9", blank_list=surgeon_level_match_trad,
                                        era="trad", sum_npg_vars= surgeon_npg_sums)

mod_match <- new_exp_surg_match_counts(data=as.data.frame(unique_surgeons_mod), hospid="newID_6", hospitals=hospitals_mod, surg_id_q ="newID_9",  blank_list=surgeon_level_match_mod,
                                       era="mod", sum_npg_vars = surgeon_npg_sums)

trad_match1 <- rcbsubset(distance.structure=trad_match, penalty=3, tol=1, exclude.penalty=NULL)
mod_match1 <- rcbsubset(distance.structure=mod_match, penalty=3, tol=1, exclude.penalty=NULL)

#### Assign Pair IDs. ####
matched_pairs_trad <- as.data.frame(trad_match1$matches)
matched_pairs_mod <- as.data.frame(mod_match1$matches)
matched_pairs_trad$case <- NA
matched_pairs_mod$case <- NA


for(i in 1:nrow(matched_pairs_trad)){
  
  temp_rowname <- as.numeric(rownames(matched_pairs_trad[i,]))
  
  matched_pairs_trad[i,c("case")] <- names(trad_match[temp_rowname])
  
  
}

for(i in 1:nrow(matched_pairs_mod)){
  
  temp_rowname <- as.numeric(rownames(matched_pairs_mod[i,]))
  
  matched_pairs_mod[i,c("case")] <- names(mod_match[temp_rowname])
  
  
}

colnames(matched_pairs_trad)[1] <- c("control")
matched_pairs_trad <- as.data.frame(matched_pairs_trad)
matched_pairs_trad$case <- as.numeric(matched_pairs_trad$case)
colnames(matched_pairs_mod)[1] <- c("control")
matched_pairs_mod <- as.data.frame(matched_pairs_mod)
matched_pairs_mod$case <- as.numeric(matched_pairs_mod$case)

trad_match_patients <- subset(as.data.frame(trad_new_exp_surg), newID_9 %in% (as.integer(matched_pairs_trad[,1])) | newID_9 %in% as.integer(matched_pairs_trad[,2]))
trad_match_patients <- as.data.frame(trad_match_patients)

mod_match_patients <- subset(as.data.frame(mod_new_exp_surg), newID_9 %in% (as.integer(matched_pairs_mod[,1])) | newID_9 %in% as.integer(matched_pairs_mod[,2]))
mod_match_patients <- as.data.frame(mod_match_patients)

patients_not_matched_trad <- subset(as.data.frame(trad_new_exp_surg), !(newID_9 %in% (as.integer(matched_pairs_trad[,1]))) & !(newID_9 %in% as.integer(matched_pairs_trad[,2])))
patients_not_matched_trad <- as.data.frame(patients_not_matched_trad)

patients_not_matched_mod <- subset(as.data.frame(mod_new_exp_surg), !(newID_9 %in% (as.integer(matched_pairs_mod[,1]))) & !(newID_9 %in% as.integer(matched_pairs_mod[,2])))
patients_not_matched_mod <- as.data.frame(patients_not_matched_mod)



trad_match_patients$pair_id_surg <- NA
for(i in 1:nrow(matched_pairs_trad)){
  
  current_new_surg <- matched_pairs_trad[i,1]
  current_exp_surg <- matched_pairs_trad[i,2]
  current_pair <- c(current_new_surg, current_exp_surg)
  
  trad_match_patients$pair_id_surg <- ifelse(trad_match_patients$newID_9 %in% current_pair, paste(i,"o",sep="_"), trad_match_patients$pair_id_surg)
  
  print(i)
  
}

mod_match_patients$pair_id_surg <- NA
for(i in 1:nrow(matched_pairs_mod)){
  
  current_new_surg <- matched_pairs_mod[i,1]
  current_exp_surg <- matched_pairs_mod[i,2]
  current_pair <- c(current_new_surg, current_exp_surg)
  
  mod_match_patients$pair_id_surg <- ifelse(mod_match_patients$newID_9 %in% current_pair, paste(i,"o",sep="_"), mod_match_patients$pair_id_surg)
  
  print(i)
  
}

##############################################
#### Create the patient-level risk score. ####

#### Create a dataset of all patients that are not eligible to be matched after the surgeon match. ####

#### Get variable names of variables that are in patients_not_matched, and add NA to the _risk_patients dataset. ####

var_names_trad <- setdiff(names(patients_not_matched_trad), names(trad_risk_patients))
var_names_mod <- setdiff(names(patients_not_matched_mod), names(mod_risk_patients))

if(length(var_names_trad)>0){
  for(i in 1:length(var_names_trad)){
    
    current_var <- var_names_trad[i]
    current_dim <- dim(trad_risk_patients)[2]
    trad_risk_patients[,current_dim+1] <- NA
    colnames(trad_risk_patients)[current_dim+1] <- current_var
    
    print(i)
    
  }
}

if(length(var_names_mod)>0){
  for(i in 1:length(var_names_mod)){
    
    current_var <- var_names_mod[i]
    current_dim <- dim(mod_risk_patients)[2]
    mod_risk_patients[,current_dim+1] <- NA
    colnames(mod_risk_patients)[current_dim+1] <- current_var
    
    print(i)
    
  }
}


risk_score_patients_trad <- rbind(trad_risk_patients, patients_not_matched_trad)
risk_score_patients_mod <- rbind(mod_risk_patients, patients_not_matched_mod)
write.csv(x=risk_score_patients_trad[,grep(x=names(risk_score_patients_trad),pattern="newID_",value=T)], file="/Users/sharpe/Desktop/Young Surgeons/data/trad_ortho_risk_patients_11.8.2018_exact_procs.csv")
write.csv(x=risk_score_patients_mod[,grep(x=names(risk_score_patients_mod),pattern="newID_",value=T)], file="/Users/sharpe/Desktop/Young Surgeons/data/mod_ortho_risk_patients_11.8.2018_exact_procs.csv")

comorbidities <- c('silb_AbdCancer',  'silb_ANGINA', 'silb_ASTHMA', 'silb_CANCER', 'silb_CHF', 'silb_CHRNLUNG', 'silb_ChrPepticUlcer',
                   'silb_COAG', 'silb_CollagenVascular', 'silb_CongenitalCoag', 'silb_CUSHINGS', 'silb_DEMENT', 'silb_DIABETES', 'silb_GRAVES', 'silb_HTN', 'silb_HYPOTHY', 'silb_LIVER',
                   'silb_PARAPLEG', 'silb_PastARR', 'silb_PastMi', 'silb_PostPulmFibr', 'silb_RENLDYS', 'silb_renlfail',
                   'silb_SEIZUR', 'silb_SickleCell', 'silb_STROKE', 'silb_THROMB', 'silb_UnstAngina', 'silb_Valvulardis', 'silb_AIDS')
aids_position <- match("silb_AIDS", comorbidities)
#### Remove AIDS from the comorbidities, since it is pretty rare.
comorbidities_for_prop_score <- comorbidities[-aids_position]



risk_vars <- c('sex_female', 'emergent_type', 'transfer_24hrs_ip', 'past_adm_6mo', ortho_npgs_for_prop_score,
               'spx_major_narrow', comorbidities_for_prop_score, 'age_70_74', 'age_75_79', 'age_80_85', 'age_85_plus', 'dead_30') 

risk_model_trad <- glm(formula = dead_30 ~ ., data=risk_score_patients_trad[,c(risk_vars)], family=binomial)

risk_model_mod_post_2010 <- glm(formula = dead_30 ~.-silb_SickleCell, data=risk_score_patients_mod[which(risk_score_patients_mod$oct_2010_plus==1),c(risk_vars[risk_vars!="silb_ChrPepticUlcer"])], family=binomial)
risk_model_mod_pre_2010 <- glm(formula = dead_30 ~.-silb_SickleCell, data=risk_score_patients_mod[which(risk_score_patients_mod$oct_2010_plus==0),c(risk_vars[risk_vars!="silb_ChrPepticUlcer"])], family=binomial)



trad_match_patients$risk_prob_trad <- predict(object=risk_model_trad, newdata=trad_match_patients, type="response")
trad_match_patients$risk_logit_trad <- predict(object=risk_model_trad, newdata=trad_match_patients, type="link", link="logit")

mod_match_patients$risk_prob_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=risk_model_mod_post_2010, newdata=mod_match_patients, type="response"),0)
mod_match_patients$risk_logit_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=risk_model_mod_post_2010, newdata=mod_match_patients, type="link", link="logit"), 0)

mod_match_patients$risk_prob_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=risk_model_mod_pre_2010, newdata=mod_match_patients, type="response"),0)
mod_match_patients$risk_logit_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=risk_model_mod_pre_2010, newdata=mod_match_patients, type="link", link="logit"), 0)

#####################################
#### Create a comorbidity score. ####

crisk_model_trad <- glm(formula = dead_30 ~., data=risk_score_patients_trad[,c(comorbidities_for_prop_score,"dead_30")], family=binomial)
#### Remove silb_ChrPepticUlcer from modern model since it produces a large estimate and standard error. ####
crisk_model_mod_post_2010 <- glm(formula=dead_30~.-silb_ChrPepticUlcer-silb_SickleCell, data=risk_score_patients_mod[which(risk_score_patients_mod$oct_2010_plus==1),c(comorbidities_for_prop_score,"dead_30")], family=binomial)
crisk_model_mod_pre_2010 <- glm(formula=dead_30~.-silb_SickleCell, data=risk_score_patients_mod[which(risk_score_patients_mod$oct_2010_plus==0),c(comorbidities_for_prop_score,"dead_30")], family=binomial)



trad_match_patients$crisk_prob_trad <- predict(object=crisk_model_trad, newdata=trad_match_patients, type="response")
trad_match_patients$crisk_logit_trad <- predict(object=crisk_model_trad, newdata=trad_match_patients, type="link", link="logit")

mod_match_patients$crisk_prob_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=crisk_model_mod_post_2010, newdata=mod_match_patients, type="response"),0)
mod_match_patients$crisk_logit_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=crisk_model_mod_post_2010, newdata=mod_match_patients, type="link", link="logit"), 0)

mod_match_patients$crisk_prob_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=crisk_model_mod_pre_2010, newdata=mod_match_patients, type="response"),0)
mod_match_patients$crisk_logit_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=crisk_model_mod_pre_2010, newdata=mod_match_patients, type="link", link="logit"), 0)



###################################
#### Create model information. ####

risk_cstat_trad <- auc(risk_score_patients_trad[,"dead_30"], risk_model_trad$fitted.values)
risk_cstat_mod_post_2010 <- auc(risk_score_patients_mod[which(risk_score_patients_mod$oct_2010_plus==1),"dead_30"], risk_model_mod_post_2010$fitted.values)
risk_cstat_mod_pre_2010 <- auc(risk_score_patients_mod[which(risk_score_patients_mod$oct_2010_plus==0),"dead_30"], risk_model_mod_pre_2010$fitted.values)

crisk_cstat_trad <- auc(risk_score_patients_trad[,"dead_30"], crisk_model_trad$fitted.values)
crisk_cstat_mod_post_2010 <- auc(risk_score_patients_mod[which(risk_score_patients_mod$oct_2010_plus==1),"dead_30"], crisk_model_mod_post_2010$fitted.values)
crisk_cstat_mod_pre_2010 <- auc(risk_score_patients_mod[which(risk_score_patients_mod$oct_2010_plus==0),"dead_30"], crisk_model_mod_pre_2010$fitted.values)

risk_model_N <- c(dim(risk_score_patients_trad)[1], dim(risk_score_patients_trad[which(risk_score_patients_trad$dead_30==1),])[1],
                  dim(risk_score_patients_mod)[1], dim(risk_score_patients_mod[which(risk_score_patients_mod$dead_30==1),])[1],
                  dim(risk_score_patients_mod[which(risk_score_patients_mod$oct_2010_plus==1),])[1], dim(risk_score_patients_mod[which(risk_score_patients_mod$oct_2010_plus==1 & risk_score_patients_mod$dead_30==1),])[1],
                  dim(risk_score_patients_mod[which(risk_score_patients_mod$oct_2010_plus==0),])[1], dim(risk_score_patients_mod[which(risk_score_patients_mod$oct_2010_plus==0 & risk_score_patients_mod$dead_30==1),])[1],
                  round(risk_cstat_trad,3), round(risk_cstat_mod_post_2010,3), round(risk_cstat_mod_pre_2010,3),
                  round(crisk_cstat_trad,3), round(crisk_cstat_mod_post_2010,3), round(crisk_cstat_mod_pre_2010,3))

names(risk_model_N) <- names(risk_model_N) <- c("Trad N","Trad N dead", "Mod N", "Mod N dead", "Post 2010 N", "Post 2010 N dead", "Pre 2010 N", "Pre 2010 N", "C-Stat Trad Risk", "C-Stat Mod Post 2010 Risk", "C-Stat Mod Pre 2010 Risk",
                                                "C-Stat Trad Comorbidity", "C-Stat Mod Post 2010 Comorbidity", "C-Stat Mod Pre 2010 Comorbidity")
risk_model_N

write.csv(x=t(risk_model_N), file="/Users/sharpe/Desktop/Young Surgeons/table/risk_model_N_ortho_11.8.2018_exact_procs.csv", row.names=T)

#################################################################
#### Create a propensity score for the patient match. ###########



prop_score_pat_vars_trad <- c('sex_female', 'emergent_type', 'transfer_24hrs_ip', 'past_adm_6mo', ortho_npgs_for_prop_score,
                              'spx_major_narrow', comorbidities_for_prop_score,'age_70_74','age_75_79','age_80_85','age_85_plus', 'new_surgeon_flag')

prop_score_pat_vars_mod <- c('sex_female', 'oct_2010_plus', 'emergent_type', 'transfer_24hrs_ip', 'past_adm_6mo', ortho_npgs_for_prop_score,
                             'spx_major_narrow', comorbidities_for_prop_score,'age_70_74','age_75_79','age_80_85','age_85_plus', 'new_surgeon_flag')





prop_score_pat_trad <- glm(formula = new_surgeon_flag ~ ., data=trad_match_patients[,c(prop_score_pat_vars_trad)], family=binomial)


prop_score_pat_mod_post_2010 <- glm(formula = new_surgeon_flag ~ .-oct_2010_plus, data=mod_match_patients[which(mod_match_patients$oct_2010_plus==1),c(prop_score_pat_vars_mod)], family=binomial)
prop_score_pat_mod_pre_2010 <- glm(formula = new_surgeon_flag ~ .-silb_CUSHINGS-oct_2010_plus-silb_SickleCell, data=mod_match_patients[which(mod_match_patients$oct_2010_plus==0),c(prop_score_pat_vars_mod)], family=binomial)

summary(prop_score_pat_trad)
summary(prop_score_pat_mod_post_2010)
summary(prop_score_pat_mod_pre_2010)


trad_match_patients$prop_prob_pat_trad <- predict(object=prop_score_pat_trad, newdata=trad_match_patients, type="response")
trad_match_patients$prop_log_pat_trad <- predict(object=prop_score_pat_trad, newdata=trad_match_patients, type="link", link="logit")

mod_match_patients$prop_prob_pat_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=mod_match_patients, type="response"), 0)
mod_match_patients$prop_log_pat_mod_post_2010 <- ifelse(mod_match_patients$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=mod_match_patients, type="link", link="logit"),0)

mod_match_patients$prop_prob_pat_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=mod_match_patients, type="response"), 0)
mod_match_patients$prop_log_pat_mod_pre_2010 <- ifelse(mod_match_patients$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=mod_match_patients, type="link", link="logit"),0)

#### Calculate the standard deviation for the propensity scores. ####

prop_score_sd_trad <- sd(trad_match_patients[, c("prop_log_pat_trad")])
prop_score_sd_mod_post_2010 <- sd(mod_match_patients[which(mod_match_patients$oct_2010_plus==1), c("prop_log_pat_mod_post_2010")])
prop_score_sd_mod_pre_2010 <- sd(mod_match_patients[which(mod_match_patients$oct_2010_plus==0), c("prop_log_pat_mod_pre_2010")])


#############################
#### Match the patients. ####


trad_match_patients <- trad_match_patients[order(-c(trad_match_patients$treatment)),]
mod_match_patients <- mod_match_patients[order(-c(mod_match_patients$treatment)),]

trad_treat <- trad_match_patients$treatment
mod_treat <- mod_match_patients$treatment


maha_vars_pat_unimp <- c(ortho_npgs, "transfer_24hrs_ip", "emergent_type", "past_adm_6mo", "silb_CHF", "silb_PastMi", "silb_PastARR", "silb_ANGINA", "silb_Valvulardis", "silb_UnstAngina",  "age_70_74", "age_75_79",
                         "age_80_85", "age_85_plus", "sex_female", "silb_RENLDYS", "silb_renlfail", "silb_LIVER", "silb_STROKE", "silb_DEMENT", "silb_SEIZUR", "silb_PARAPLEG", "silb_CANCER", "silb_AbdCancer", "silb_CHRNLUNG",
                         "silb_ASTHMA", "silb_PostPulmFibr", "silb_DIABETES", "silb_HTN", "jan_2012_plus",  "spx_major_narrow", "silb_CongenitalCoag", "silb_THROMB", "silb_CollagenVascular",
                         "silb_AIDS", "silb_COAG", "silb_SickleCell", "silb_HYPOTHY", "silb_ChrPepticUlcer", "silb_CUSHINGS", "silb_GRAVES", "admsn_in_quarters")


maha_vars_pat_imp <- c(ortho_npgs, "transfer_24hrs_ip", "emergent_type", "past_adm_6mo", "silb_CHF", "silb_PastMi", "silb_PastARR", "silb_ANGINA", "silb_Valvulardis", "silb_UnstAngina",  "age_70_74", "age_75_79",
                       "age_80_85", "age_85_plus", "sex_female", "silb_RENLDYS", "silb_renlfail", "silb_LIVER", "silb_STROKE", "silb_DEMENT", "silb_SEIZUR", "silb_PARAPLEG", "silb_CANCER", "silb_AbdCancer", "silb_CHRNLUNG",
                       "silb_ASTHMA", "silb_PostPulmFibr", "silb_DIABETES", "admsn_in_quarters")

trad_match_patients$exact_group_2 <- paste(trad_match_patients$pair_id_surg, trad_match_patients$Narrow_Procedure_Group, sep="_")
mod_match_patients$exact_group_2 <- paste(mod_match_patients$pair_id_surg, mod_match_patients$oct_2010_plus, mod_match_patients$Narrow_Procedure_Group, sep="_")



exact_vars_pat <- c("exact_group_2")
caliper_var_pat_trad <- c("prop_log_pat_trad","risk_logit_trad","spx_major_narrow", "silb_CHRNLUNG", "silb_DIABETES", "comorb_sum", "crisk_logit_trad",
                          "silb_CollagenVascular", "age_85_plus", "emergent_type", "transfer_24hrs_ip", "silb_PARAPLEG", "silb_CANCER", "silb_ASTHMA", "parapleg_mod", "congcoag_mod",
                          "collvasc_mod", "collvasc_trad", "hypo_trad", "silb_ANGINA",
                          "silb_PastMi","silb_UnstAngina","silb_STROKE","silb_SEIZUR","silb_renlfail", "renal_fail_trad", "hypo_mod","unstang_mod","seizure_mod","parapleg_mod","stroke_trad","past_adm_trad")

caliper_var_pat_mod <- c("prop_log_pat_mod_post_2010","prop_log_pat_mod_pre_2010","risk_logit_mod_post_2010","risk_logit_mod_pre_2010","spx_major_narrow", "silb_CHRNLUNG", "silb_DIABETES", "comorb_sum",
                         "crisk_logit_mod_post_2010","crisk_logit_mod_pre_2010", "silb_CollagenVascular", "age_85_plus", "emergent_type", "transfer_24hrs_ip", "silb_PARAPLEG", "silb_CANCER", "silb_ASTHMA",
                         "parapleg_mod", "congcoag_mod", "collvasc_mod", "collvasc_trad", "hypo_trad", "silb_ANGINA",
                         "silb_PastMi","silb_UnstAngina","silb_STROKE","silb_SEIZUR","silb_renlfail", "renal_fail_trad", "hypo_mod","unstang_mod","seizure_mod","parapleg_mod","stroke_trad","past_adm_trad")

near_exact_vars <- NA



mnf.names <- vector("list",length=3)
# mnf.names[[1]] <- c(ortho_npgs)
mnf.names[[1]] <- c("emergent_type", "transfer_24hrs_ip")
mnf.names[[2]] <- c(mnf.names[[1]],  "silb_CHF", "silb_PastMi", "silb_PastARR", "silb_ANGINA", "silb_Valvulardis", "silb_UnstAngina", "silb_HTN", "silb_RENLDYS","silb_renlfail","silb_LIVER","silb_STROKE", "silb_DEMENT", "silb_SEIZUR", "silb_PARAPLEG","silb_CHRNLUNG")
mnf.names[[3]] <- c(mnf.names[[2]], "age_65_69", "age_70_74", "age_75_79", "age_80_85", "age_85_plus","sex_female", "past_adm_6mo", "jan_2012_plus")


maha_data_pat_unimp_trad <- trad_match_patients[,maha_vars_pat_unimp]
maha_data_pat_unimp_mod <- mod_match_patients[,maha_vars_pat_unimp]

maha_data_pat_imp_trad <- trad_match_patients[,maha_vars_pat_imp]
maha_data_pat_imp_mod <- mod_match_patients[,maha_vars_pat_imp]
exact_matrix_pat_trad <- trad_match_patients[,exact_vars_pat]
exact_matrix_pat_mod <- mod_match_patients[,exact_vars_pat]
calip_data_pat_trad <- trad_match_patients[,caliper_var_pat_trad]
calip_data_pat_mod <- mod_match_patients[,caliper_var_pat_mod]



#### Create the unimportant distance matrices, divide by its mean and then add penalties. ####
dist_struc_pats_unimp_trad <- build.dist.struct(z=trad_treat, X = maha_data_pat_unimp_trad, exact=exact_matrix_pat_trad, calip.option="none")

avg_unimp_trad <- mean(unlist(dist_struc_pats_unimp_trad))

for(i in 1:length(dist_struc_pats_unimp_trad)){
  for(j in 1:length(dist_struc_pats_unimp_trad[[i]])){
    dist_struc_pats_unimp_trad[[i]][j] <- (dist_struc_pats_unimp_trad[[i]][j]/avg_unimp_trad)*1000
  }
  names(dist_struc_pats_unimp_trad[[i]]) <- names(dist_struc_pats_unimp_trad[[i]])
}

####

dist_struc_pats_unimp_mod <- build.dist.struct(z=mod_treat, X = maha_data_pat_unimp_mod, exact=exact_matrix_pat_mod, calip.option="none")

avg_unimp_mod <- mean(unlist(dist_struc_pats_unimp_mod))

for(i in 1:length(dist_struc_pats_unimp_mod)){
  for(j in 1:length(dist_struc_pats_unimp_mod[[i]])){
    dist_struc_pats_unimp_mod[[i]][j] <- (dist_struc_pats_unimp_mod[[i]][j]/avg_unimp_mod)*1000
  }
  names(dist_struc_pats_unimp_mod[[i]]) <- names(dist_struc_pats_unimp_mod[[i]])
}

####

dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=maha_data_pat_unimp_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", pen.size=2000, var="admsn_in_quarters")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc  =dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_trad,cal.size=0.20, pen.size=1000, var="prop_log_pat_trad")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=5500, var="risk_logit_trad")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=5500, var="crisk_logit_trad")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=2500, var="silb_DIABETES")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=2500, var="silb_CollagenVascular")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="control", pen.size=2500, var="silb_CHRNLUNG")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=2500, var="silb_ASTHMA")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=5500, var="comorb_sum")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=1000, var="hypo_trad")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=7000, var="renal_fail_trad")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=4000, var="past_adm_trad")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=4000, var="stroke_trad")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_CHRNLUNG")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_PARAPLEG")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_DIABETES")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_CANCER")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_ASTHMA")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_CollagenVascular")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", pen.size=2000, var="comorb_sum")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="age_85_plus")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=5000, var="emergent_type")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=5000, var="transfer_24hrs_ip")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="spx_major_narrow")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="collvasc_trad")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_ANGINA")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_PastMi")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_UnstAngina")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_STROKE")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_SEIZUR")
dist_struc_pats_unimp_trad <- apply.penalty(dist.struc = dist_struc_pats_unimp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_renlfail")

####

dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=maha_data_pat_unimp_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", pen.size=2000, var="admsn_in_quarters")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc  =dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_mod_post_2010,cal.size=0.20, pen.size=1000, var="prop_log_pat_mod_post_2010")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc  =dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_mod_pre_2010,cal.size=0.20, pen.size=1000, var="prop_log_pat_mod_pre_2010")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=5500, var="risk_logit_mod_post_2010")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=5500, var="risk_logit_mod_pre_2010")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=5500, var="crisk_logit_mod_post_2010")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=5500, var="crisk_logit_mod_pre_2010")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=2500, var="silb_DIABETES")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=2500, var="silb_CollagenVascular")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="control", pen.size=2500, var="silb_CHRNLUNG")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=2500, var="silb_ASTHMA")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=5500, var="comorb_sum")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=4000, var="parapleg_mod")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=4000, var="seizure_mod")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=4000, var="unstang_mod")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="control", pen.size=4000, var="hypo_mod")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_CHRNLUNG")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_PARAPLEG")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_DIABETES")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_CANCER")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_ASTHMA")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_CollagenVascular")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", pen.size=2000, var="comorb_sum")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="age_85_plus")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=5000, var="emergent_type")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=5000, var="transfer_24hrs_ip")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="spx_major_narrow")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="parapleg_mod")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="congcoag_mod")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="collvasc_mod")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_ANGINA")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_PastMi")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_UnstAngina")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_STROKE")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_SEIZUR")
dist_struc_pats_unimp_mod <- apply.penalty(dist.struc = dist_struc_pats_unimp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_renlfail")


################################################################
#### Create the important distance matrices; add penalties. ####

dist_struc_pats_imp_trad <- build.dist.struct(z=trad_treat, X = maha_data_pat_imp_trad, exact=exact_matrix_pat_trad, calip.option="none")

avg_imp_trad <- mean(unlist(dist_struc_pats_imp_trad))

for(i in 1:length(dist_struc_pats_imp_trad)){
  for(j in 1:length(dist_struc_pats_imp_trad[[i]])){
    dist_struc_pats_imp_trad[[i]][j] <- (dist_struc_pats_imp_trad[[i]][j]/avg_imp_trad)*1500
  }
  names(dist_struc_pats_imp_trad[[i]]) <- names(dist_struc_pats_imp_trad[[i]])
}

max(unlist(dist_struc_pats_imp_trad))

####

dist_struc_pats_imp_mod <- build.dist.struct(z=mod_treat, X = maha_data_pat_imp_mod, exact=exact_matrix_pat_mod, calip.option="none")

avg_imp_mod <- mean(unlist(dist_struc_pats_imp_mod))

for(i in 1:length(dist_struc_pats_imp_mod)){
  for(j in 1:length(dist_struc_pats_imp_mod[[i]])){
    dist_struc_pats_imp_mod[[i]][j] <- (dist_struc_pats_imp_mod[[i]][j]/avg_imp_mod)*1500
  }
  names(dist_struc_pats_imp_mod[[i]]) <- names(dist_struc_pats_imp_mod[[i]])
}

max(unlist(dist_struc_pats_imp_mod))


####

dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=maha_data_pat_imp_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.5,  pen.size=2000, var="admsn_in_quarters")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc  =dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_trad, cal.size=0.20, pen.size=2000, var="prop_log_pat_trad")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc  =dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=4500, var="risk_logit_trad")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc  =dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=4500, var="crisk_logit_trad")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc  =dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=3500, var="silb_DIABETES")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc  =dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=3500, var="silb_CollagenVascular")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc  =dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="control", pen.size=3500, var="silb_CHRNLUNG")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc  =dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=3500, var="silb_ASTHMA")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc  =dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=5000, var="comorb_sum")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc  =dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=2000, var="hypo_trad")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc  =dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=7000, var="renal_fail_trad")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=3000, var="silb_CHRNLUNG")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="silb_PARAPLEG")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="silb_DIABETES")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="silb_CANCER")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="silb_CollagenVascular")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", pen.size=2000, var="comorb_sum")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=2000, var="age_85_plus")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="emergent_type")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="transfer_24hrs_ip")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="collvasc_trad")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_ANGINA")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_PastMi")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_UnstAngina")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_STROKE")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_SEIZUR")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_renlfail")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=4000, var="past_adm_trad")
dist_struc_pats_imp_trad <- apply.penalty(dist.struc = dist_struc_pats_imp_trad, X=calip_data_pat_trad, exact=exact_matrix_pat_trad, z=trad_treat, type="directional",direction="case", pen.size=4000, var="stroke_trad")

####

dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=maha_data_pat_imp_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.5,  pen.size=2000, var="admsn_in_quarters")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_mod_post_2010, cal.size=0.20, pen.size=2000, var="prop_log_pat_mod_post_2010")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=TRUE, own.sd=TRUE, sd.var=prop_score_sd_mod_pre_2010, cal.size=0.20, pen.size=2000, var="prop_log_pat_mod_pre_2010")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=4500, var="risk_logit_mod_post_2010")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=4500, var="risk_logit_mod_pre_2010")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=4500, var="crisk_logit_mod_post_2010")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=4500, var="crisk_logit_mod_pre_2010")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=3500, var="silb_DIABETES")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=3500, var="silb_CollagenVascular")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="control", pen.size=3500, var="silb_CHRNLUNG")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=3500, var="silb_ASTHMA")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc  =dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=5000, var="comorb_sum")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=3000, var="silb_CHRNLUNG")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="silb_PARAPLEG")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="silb_DIABETES")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="silb_CANCER")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="silb_CollagenVascular")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", pen.size=2000, var="comorb_sum")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=2000, var="age_85_plus")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="emergent_type")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="transfer_24hrs_ip")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="parapleg_mod")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="congcoag_mod")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="caliper", cal.sd=FALSE, cal.size=0.5,  pen.size=5000, var="collvasc_mod")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_ANGINA")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_PastMi")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_UnstAngina")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=2000, var="silb_STROKE")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_SEIZUR")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="scaled", cal.sd=FALSE, cal.size=0.50, pen.size=1000, var="silb_renlfail")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=4000, var="parapleg_mod")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=4000, var="seizure_mod")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="case", pen.size=4000, var="unstang_mod")
dist_struc_pats_imp_mod <- apply.penalty(dist.struc = dist_struc_pats_imp_mod, X=calip_data_pat_mod, exact=exact_matrix_pat_mod, z=mod_treat, type="directional",direction="control", pen.size=4000, var="hypo_mod")

####



dist_struc_pats_trad <- vector("list", length(dist_struc_pats_unimp_trad))

for(i in 1:length(dist_struc_pats_unimp_trad)){
  
  
  for(j in 1:length(dist_struc_pats_unimp_trad[[i]])){
    
    dist_struc_pats_trad[[i]][j] <- (dist_struc_pats_unimp_trad[[i]][j]) + (dist_struc_pats_imp_trad[[i]][j]) 
    
  }
  
  names(dist_struc_pats_trad[[i]]) <- names(dist_struc_pats_unimp_trad[[i]])
}


####

dist_struc_pats_mod <- vector("list", length(dist_struc_pats_unimp_mod))

for(i in 1:length(dist_struc_pats_unimp_mod)){
  
  
  for(j in 1:length(dist_struc_pats_unimp_mod[[i]])){
    
    dist_struc_pats_mod[[i]][j] <- (dist_struc_pats_unimp_mod[[i]][j]) + (dist_struc_pats_imp_mod[[i]][j]) 
    
  }
  
  names(dist_struc_pats_mod[[i]]) <- names(dist_struc_pats_unimp_mod[[i]])
}

####

treated_info_pats_trad <- subset(trad_match_patients, treatment==1)
control_info_pats_trad <- subset(trad_match_patients, treatment==0)

treated_info_pats_mod <- subset(mod_match_patients, treatment==1)
control_info_pats_mod <- subset(mod_match_patients, treatment==0)


patient_match_trad <- rcbsubset(distance.structure = dist_struc_pats_trad, fb.list = mnf.names, treated.info = treated_info_pats_trad, control.info = control_info_pats_trad,
                                penalty=2, tol=10)

patient_match_mod <- rcbsubset(distance.structure = dist_struc_pats_mod, fb.list = mnf.names, treated.info = treated_info_pats_mod, control.info = control_info_pats_mod,
                               penalty=2, near.exact=c(ortho_npgs), tol=10)


tem_match_pats_trad <- treated_info_pats_trad[as.numeric(rownames(patient_match_trad$matches)),]
sda_match_pats_trad <- control_info_pats_trad[(patient_match_trad$matches),]
pair_id_pats_trad <- c(seq(1:dim(tem_match_pats_trad)[1]), seq(1:dim(sda_match_pats_trad)[1]))
combined_pats_trad <- rbind(tem_match_pats_trad, sda_match_pats_trad)
combined_pats_trad$pair_id_pats <- pair_id_pats_trad

tem_match_pats_mod <- treated_info_pats_mod[as.numeric(rownames(patient_match_mod$matches)),]
sda_match_pats_mod <- control_info_pats_mod[(patient_match_mod$matches),]
pair_id_pats_mod <- c(seq(1:dim(tem_match_pats_mod)[1]), seq(1:dim(sda_match_pats_mod)[1]))
combined_pats_mod <- rbind(tem_match_pats_mod, sda_match_pats_mod)
combined_pats_mod$pair_id_pats <- pair_id_pats_mod

#### Add propensity scores to the before datasets.

NEW_PATIENTS_BEFORE$prop_prob_pat_trad <- ifelse(NEW_PATIENTS_BEFORE$era_trad==1, predict(object=prop_score_pat_trad, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_trad <- ifelse(NEW_PATIENTS_BEFORE$era_trad==1, predict(object=prop_score_pat_trad, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)
NEW_PATIENTS_BEFORE$prop_prob_pat_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)
NEW_PATIENTS_BEFORE$prop_prob_pat_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$prop_log_pat_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)



exp_patients_before$prop_prob_pat_trad <- ifelse(exp_patients_before$era_trad==1, predict(object=prop_score_pat_trad, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_trad <- ifelse(exp_patients_before$era_trad==1, predict(object=prop_score_pat_trad, newdata=exp_patients_before, type="link", link="logit"),0)
exp_patients_before$prop_prob_pat_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=prop_score_pat_mod_post_2010, newdata=exp_patients_before, type="link", link="logit"),0)
exp_patients_before$prop_prob_pat_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$prop_log_pat_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=prop_score_pat_mod_pre_2010, newdata=exp_patients_before, type="link", link="logit"),0)

#### Add risk scores to the before datasets.

NEW_PATIENTS_BEFORE$risk_prob_trad <- ifelse(NEW_PATIENTS_BEFORE$era_trad==1, predict(object=risk_model_trad, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$risk_logit_trad <- ifelse(NEW_PATIENTS_BEFORE$era_trad==1,predict(object=risk_model_trad, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)
NEW_PATIENTS_BEFORE$risk_prob_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1,predict(object=risk_model_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$risk_logit_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1,predict(object=risk_model_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)
NEW_PATIENTS_BEFORE$risk_prob_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0,predict(object=risk_model_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$risk_logit_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0,predict(object=risk_model_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)


exp_patients_before$risk_prob_trad <- ifelse(exp_patients_before$era_trad==1,predict(object=risk_model_trad, newdata=exp_patients_before, type="response"),0)
exp_patients_before$risk_logit_trad <- ifelse(exp_patients_before$era_trad==1,predict(object=risk_model_trad, newdata=exp_patients_before, type="link", link="logit"),0)
exp_patients_before$risk_prob_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=risk_model_mod_post_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$risk_logit_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=risk_model_mod_post_2010, newdata=exp_patients_before, type="link", link="logit"),0)
exp_patients_before$risk_prob_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=risk_model_mod_pre_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$risk_logit_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=risk_model_mod_pre_2010, newdata=exp_patients_before, type="link", link="logit"),0)


#### Add comorbidity risk scores to the before datasets. ####

NEW_PATIENTS_BEFORE$crisk_prob_trad <- ifelse(NEW_PATIENTS_BEFORE$era_trad==1, predict(object=crisk_model_trad, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$crisk_logit_trad <- ifelse(NEW_PATIENTS_BEFORE$era_trad==1,predict(object=crisk_model_trad, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)
NEW_PATIENTS_BEFORE$crisk_prob_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1,predict(object=crisk_model_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$crisk_logit_mod_post_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==1,predict(object=crisk_model_mod_post_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)
NEW_PATIENTS_BEFORE$crisk_prob_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0,predict(object=crisk_model_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="response"),0)
NEW_PATIENTS_BEFORE$crisk_logit_mod_pre_2010 <- ifelse(NEW_PATIENTS_BEFORE$era_mod==1 & NEW_PATIENTS_BEFORE$oct_2010_plus==0,predict(object=crisk_model_mod_pre_2010, newdata=NEW_PATIENTS_BEFORE, type="link", link="logit"),0)

exp_patients_before$crisk_prob_trad <- ifelse(exp_patients_before$era_trad==1,predict(object=crisk_model_trad, newdata=exp_patients_before, type="response"),0)
exp_patients_before$crisk_logit_trad <- ifelse(exp_patients_before$era_trad==1,predict(object=crisk_model_trad, newdata=exp_patients_before, type="link", link="logit"),0)
exp_patients_before$crisk_prob_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=crisk_model_mod_post_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$crisk_logit_mod_post_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==1, predict(object=crisk_model_mod_post_2010, newdata=exp_patients_before, type="link", link="logit"),0)
exp_patients_before$crisk_prob_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=crisk_model_mod_pre_2010, newdata=exp_patients_before, type="response"),0)
exp_patients_before$crisk_logit_mod_pre_2010 <- ifelse(exp_patients_before$era_mod==1 & exp_patients_before$oct_2010_plus==0, predict(object=crisk_model_mod_pre_2010, newdata=exp_patients_before, type="link", link="logit"),0)


NEW_PATIENTS_BEFORE$treatment <- 1



###############################################
#### Save necessary pieces of information. ####

#### Save the before match datasets. ####
write.csv(x=NEW_PATIENTS_BEFORE, file="/Users/sharpe/Desktop/Young Surgeons/data/NEW_PATIENTS_BEFORE_ortho_11.8.2018_exact_procs.csv")
write.csv(x=exp_patients_before, file="/Users/sharpe/Desktop/Young Surgeons/data/exp_patients_before_ortho_11.8.2018_exact_procs.csv")



#### Save the matched data. ####
write.csv(combined_pats_trad,  "/Users/sharpe/Desktop/Young Surgeons/data/simul_trad_patient_match_ortho_11.8.2018_exact_procs.csv")
write.csv(combined_pats_mod,  "/Users/sharpe/Desktop/Young Surgeons/data/simul_mod_patient_match_ortho_11.8.2018_exact_procs.csv")

#### Save the models. ####

sumrisk_trad <- summary(risk_model_trad)
write.csv(sumrisk_trad$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/table/risk_model_trad_ortho_11.8.2018_exact_procs.csv")

sumrisk_mod <- summary(risk_model_mod_pre_2010)
write.csv(sumrisk_mod$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/table/risk_model_mod_ortho_pre_2010_11.8.2018_exact_procs.csv")

sumrisk_mod <- summary(risk_model_mod_post_2010)
write.csv(sumrisk_mod$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/table/risk_model_mod_ortho_post_2010_11.8.2018_exact_procs.csv")


sumcrisk_trad <- summary(crisk_model_trad)
write.csv(sumcrisk_trad$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/table/crisk_model_trad_ortho_11.8.2018_exact_procs.csv")

sumcrisk_mod <- summary(crisk_model_mod_pre_2010)
write.csv(sumcrisk_mod$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/table/crisk_model_mod_ortho_pre_2010_11.8.2018_exact_procs.csv")
sumcrisk_mod <- summary(crisk_model_mod_post_2010)
write.csv(sumcrisk_mod$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/table/crisk_model_mod_ortho_post_2010_11.8.2018_exact_procs.csv")


prop_score_trad_summary <- summary(prop_score_pat_trad)
write.csv(prop_score_trad_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/table/propensity_model_trad_ortho_11.8.2018_exact_procs.csv")

prop_score_mod_summary <- summary(prop_score_pat_mod_pre_2010)
write.csv(prop_score_mod_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/table/propensity_model_mod_ortho_pre_2010_11.8.2018_exact_procs.csv")
prop_score_mod_summary <- summary(prop_score_pat_mod_post_2010)
write.csv(prop_score_mod_summary$coefficients, file="/Users/sharpe/Desktop/Young Surgeons/table/propensity_model_mod_ortho_post_2010_11.8.2018_exact_procs.csv")








